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Abstract 

The Ebola virus is spreading throughout West Africa and is causing thousands of deaths. In 
order to quantify the effectiveness of different strategies for controlling the spread, we develop 
a mathematical model in which the propagation of the Ebola virus through Liberia is caused 
by travel between counties. For the initial months in which the Ebola virus spreads, we find 
that the arrival times of the disease into the counties predicted by our model are compatible 
with World Health Organization data, but we also find that reducing mobility is insufficient 
to contain the epidemic because it delays the arrival of Ebola virus in each county by only a 
few weeks. We study the effect of a strategy in which safe burials are increased and effective 
hospitalisation instituted under two scenarios: (i) one implemented in mid-July 2014 and (ii) 
one in mid-August—which was the actual time that strong interventions began in Liberia. We 
find that if scenario (i) had been pursued the lifetime of the epidemic would have been three 
months shorter and the total number of infected individuals 80% less than in scenario (ii). Our 
projection under scenario (ii) is that the spreading will stop by mid-spring 2015. 
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INTRODUCTION 


For a fleeting moment last spring, the epidemic sweeping West Africa might have been 
stopped. But the opportunity to control the virus, which has now caused more than 7,800 
deaths, was lost [l|. 

The current Ebola outbreak in Western Africa is one of the deadliest and most persis¬ 


tent of epidemics [2j. According to World Health Organization data |3[ as of 31 December 
2014 there have been 20,171 cases and 7,889 deaths in three countries alone: Guinea, 
Sierra Leone, and Liberia. These numbers increase when cases and deaths from countries 
in which the outbreak has been officially declared over 4] are included. 

Cultural, economic, and political factors in that region of Western Africa 
hampered the effectiveness of the intervention strategies used by the health authorities. 
Because of a lack of reliable information about local patterns of the spreading of the Ebola 


have 


virus disease (EVD) 10H12]. the strategies currently being used, including the mobilisation 


of resources, the creation of new Ebola treatment centers (ETC), the development of safe 
burial procedures, and the international coordination of the efforts jul as of 1 January 
2015 have been only partially successful. 

Legrand et al. 14!] developed a seminal mathematical stochastic model with full mixing 
that reproduces the 1995 EVD outbreak in the Congo and the 2000 outbreak in Uganda. 
The population is divided into six compartments. Individuals in the susceptible compart¬ 
ment transition to exposed compartment and to the infectious compartment when they 
become infected. A percentage of these infected individuals are hospitalised and there 
are two possible outcomes: (i) they die, but before they are removed from the epidemic 
system they transition into the funeral compartment and infect other susceptible individ¬ 
uals, or (ii) they are removed from the system because they are cured. The maximum 
likelihood method is used to calibrate the model with the data. 

Rivers et al. [15] used a deterministic version of this model and least-squares optimi¬ 


sation to fit the current Liberia and Sierra Leone outbreak data. Their model indicated 
that the epidemic would not reach its peak until 31 December 2015. Gomes et al. 116] 

n 

estimated the transmission coefficients using the model provided by Ref. [14], a Global 
Epidemic and Mobility model that uses a structured metapopulation scheme, integrating 
the stochastic modelling of the disease dynamic, high resolution census and human rno- 
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bility patterns at the global scale using a high resolution population data 17], [18]. The 
parameters were estimated by fitting the total number of cumulative deaths from Liberia, 
Sierra Leone, and Guinea during the period 6 July - 9 August 2014. The transmission 
parameters obtained were used to forecast three months of EVD propagation in West 
Africa and the probability of its spreading internationally. They found that the risk of 
cases spreading to other countries was low. Poletto et al. [3] used the same model and 
found that reducing the number of travellers crossing international boundaries delays the 


arrival of EVD by only a few weeks. Merler et al. 


Ref. 


20] used methods similar to those in 


to model the effect of epidemic spreading between geographical regions. They 
took into account the movements of non-infected individuals who were assisting in health¬ 
care facilities, those who took care of non-hospitalised infected individuals, and those who 
attended funerals. 

Population mobility—the movement of individuals seeking safer areas, better health 
infrastructures, or food supplies—strongly affects disease propagation and plays a major 
role in epidemic spreading and in the effectiveness of any intervention scheme 2l|. In 
Liberia, 54% of the population over the age of 14 are internally displaced 22], Under¬ 
standing these patterns of movement is essential when planning interventions to contain 
regional outbreaks. In recent years a number of mobility studies have been published 


23- 


251 ]. including Wesolowski et al. 2l|, who used mobile telephone network data to analyse 


mobility patterns that could be useful to understand the Ebola outbreak. They analysed 
data sources from mobile phone call records (CDRs), national census microdata samples, 
and spatial population data in order to estimate domestic and international mobility pat¬ 
terns in West African countries. The best mobility estimates were obtained for Senegal, 
Cote d’Ivoire, and Kenya, and Wesolowski et al. |2l|, l25[] used them to produce a spatial 
interaction model of national mobility patterns in order to estimate how the EVD affected 
regions are connected by population flows. 

We use a stochastic compartmental model and a set of differential equations, which are 
the quasi-deterministic representation of a stochastic model, to understand how popula¬ 
tion mobility affects the spreading of EVD between regions (comities) within Liberia. Our 
model quantifies how mobility between comities affects epidemic spreading inside Liberia, 
and we find that although reducing mobility among comities delays the spread of Ebola, 
it does not contain it. Our study indicates that the response implemented in August 2014 
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will result the extinction of the epidemic by mid-spring 2015, but it also indicates that 
an earlier response would have been extremely effective in containing the disease. 


RESULTS 


Model 


In our model we classify individuals as susceptible (S), exposed (E), i.e., infected but 
not infectious, infected (I), hospitalised (H), recovered (R), i.e., either cured or dead with 
a safe burial that does not transmit the disease, or dead (F) with an unsafe burial that 
transmits the disease. We also classify infected and hospitalised individuals according 
to their fate: those who are infected, will be hospitalised, and will die (Idh), those who 
are infected, won’t be hospitalised, and will die (Idnh), those who are infected, will be 
hospitalised, and will recover (Irh), those who are infected, won’t be hospitalised, and 
will recover (Irnh), those who are hospitalised and will die (Hd), and those who are 
hospitalised and will recover (Hr). The symbols S, E, I, H, R, and F indicate both the 
classification and the population percentage within the classification. 

Figure |T| shows a schematic presentation of the model indicating the compartmental 
states (red boxes) and the transition rates among the states (connecting arrows). The /, 
/ = Idh+Idnh+Irh+Irnh represents the total number of infected individuals, and H = 
Hr + Hr, the total number of those hospitalised. Table [T| shows the different parameters 
used to calculate the transition rates among the different compartmental states, and 
Table SI (see Supplementary Information) shows the N T = 12 transitions between states 
and their rates A* with % — 1, ...Nr. 

To determine how geographic mobility spreads the disease, we utilise the model of West 


21 


251. In their 


African regional transportation patterns developed by Wesolowski et al. 
research they applied a gravity model to mobile phone data for Senegal to estimate the flow 
of individuals between counties in Liberia. Although these movement data are “historical” 
and do not reflect how local population behaviour may have changed in response to the 
current crisis, we assume the patterns of mobility obtained in the Wesolowski model 
J2I, 2fj] still represent a good approximation of the routine commuting patterns of the 
population in Liberia prior to the outbreak. This is different to post-outbreak models 
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FIG. 1: A schematic of the transitions between different states of our model for 
the EVD spreading in West Africa 2014 and their respective transition rates. In 

the model, the population is divided into ten compartmental states (See Table SI): Susceptible 
(S) individuals who in contact with infected individuals can become exposed (E). These E 
individuals after the incubation period become infected and follow four different scenarios: (i) 
Infected individuals that will be cured —recovered— without hospitalisation (Irnh); (h) Infected 
individuals who will be cured (Irh) after spending a period on a hospital (H r); (iii) Infected 
individuals without being hospitalised (Idnh) who will die and may infect other individuals 
in their funerals (F); and (iv) Infected individuals (Idh) that even after spending a period 
in a hospital (Hp) will die and may also spread the infection in the funerals ( F ). Recovered 
individual (R), can be cured or dead. 


that describe travel patterns that reflect human efforts to avoid the disease or to attend 
funerals of epidemic victims (see Ref. 20] and references in therein). 

We assume that there is a flow of individuals between all N co = 15 counties of Liberia, 
and that only susceptible or exposed individuals can travel between counties. Thus the 
deterministic evolution equations for the number of individuals in each state in county c 
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in our model are 
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where a c is the total rate of mobility in each county c and is given by 
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where x ( x c ) is the number of individuals (susceptible or exposed), in county j (c), Nj 
(N c ) the total population of county j (c), and r 3 c and r C j the mobility rates from county 
j -A c and from county c -A j, respectively. Note that due to mobility the population 
in each county changes, but since this evolution is much slower than the dynamics of 
the disease spreading, we consider N c to be constant (in our model without restriction 
on the mobility, the population in each county changes less than a 5% each year). In 
addition, in this model we disregard the mobility inside each county, i.e., we assume that 
the population is fully mixed. Because recovered individuals are unable to transmit the 
disease or be reinfected, they do not affect the results of our model and we disregard their 
movements between counties. 

In the context of complex network research, this model of mobility between counties 
breaks the traditional full-mixing a ppr oach because each county can be thought of as a 


node of a metapopulation network 


27] in which the weight of each link is proportional 


6 













to the mobility flow. Note that if in Eqs. (TTj MIlOl) we drop the index c and disregard the 
flow mobility we are no longer taking the comities into account, and we have a scenario 
that represents the spread throughout the entire country. 


Transmission rates estimated 


01 . 


According to WHO data [3|, the first index case (patient zero) was diagnosed in Lofa on 
17 March 2014. Thus our initial conditions in Lofa are (i) one infected individual in that 
county and (ii) the rest of the population susceptible. The estimated rates of transmission 
in day^ 1 obtained (using the method presented in the section Methods: Calibration with 
the deterministic equations ^ are fdi = 0.14 [0,0.26] in the community, /3h = 0.29 [0,0.92] 
in the hospitals, and /3p = 0.40 [0, 0.99] at the funerals, where the intervals correspond 
to the values used to obtain the average rates of transmission obtained from the Akaike 


criterion. From these rates, we construct the next-generation matrix 28,|29[] (see Methods: 


Estimation of Ro) in order to compute the reproductive number Ro, defined as the average 
number of people in a susceptible population one infected individual infects during his or 
her infectious period. This parameter is fundamental when predicting whether a disease 
can reach a macroscopic fraction of individuals jsoj]. For a critical value R 0 = 1 there is 
a phase transition below which no epidemic takes place, and the disease is only a small 
outbreak, while for Rq > 1 the probability that an epidemic spreading develops is greater 


than zero 


30]. For the values of rates of transmission given above, we find that the 


reproductive number of the current EVD outbreak is i? 0 = 2.11 [1.88,2.71], well above 
the critical threshold Ro — 1, where the interval was obtained from the transmission 
coefficient selected from the Akaike criterion. This value of R 0 is compatible with the one 


obtained in Ref. m . We run our stochastic simulations presented in Methods: Stochastic 
model for these estimated values in order to compare the total number of cases with the 
data given by WHO before the interventions began in the middle of August 2014 13]. 
Figure [21(a) plots the number of cumulative cases as a function of time for 1000 realisations 


of our stochastic model and compares the results with WHO data [3[ without any shift 
correction. The individual realisations have the same shape as the data but due to the 
stochasticity at the beginning of the outbreak the exponential increase in the number of 
cases occurs at different moments. 
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FIG. 2: Cumulative number of cases in Liberia for the parameters given in Table |TJ 

Cumulative number of cases obtained with our stochastic model with the transition presented 
in Table SI and Eq. dill) in Liberia with 1000 realisations (gray lines) and the data (symbols) 
without temporal shift (a) and (b) with a temporal shift using s c = 200. The transmission 
coefficients /?/ = 0.14, (3h = 0.29 and /3p = 0.40 were obtained as explained in Methods: 
Calibration with the deterministic equations. From the WHO’s data the index case is located at 
Lofa on March 17 2014. 

Figure [21(b) plots the cumulative number of cases as a function of time with the initial 
conditions explained above when a temporal shift is applied to the stochastic simula¬ 
tions. The agreement between the simulations and the data indicates that our model 
can successfully represent the dynamics of the spreading of the current Ebola outbreak in 
Liberia. 

The geographical spread of Ebola cases across Liberia due to mobility 

The mobility among the 15 counties allows us to compute the arrival time t a in each 
county, assuming that the index case was in Lofa on 17 March 2014. Figure [3] shows the 
violin plots of the arrival times t a of the disease as it spreads from Lofa County into the 
other 14 Liberian counties and compares our results with those supplied in the WHO 
reports (circles). 
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FIG. 3: Figure a: Violin plots representing the distribution o: 
considering the mobility flow of individuals among counties 


arrival time t a to each county 


21] without any restriction on the 


mobility. The results are obtained from our stochastic model with the estimated transmission 
coefficients over 1000 realisations. Error bars indicate the 95% confidence interval. From WHO’s 
reports the index case (patient zero) was located at Lofa at 17 of March 2014. The circles 
represent the values of t a reported by WHO. The very early case of Margibi is below the 5% 


probability, and it is explained in Ref, 


'• 31]- 


Figure b: Violin plots representing the distribution 


of arrival time t a to each county reducing the mobility among counties by 80%. 


Comparing the results of our predictions of the arrival times of the first case as it 
spreads to the other counties with the WHO data (see Fig. Ek), all counties except Margibi 
and Grand Gedeh fall into a 95% confidence interval. This could be caused by (i) an 

n 

underestimation of the number of cases in the WHO data [3j due to a lack of information 
[l|, or (ii) because the data recorded are actually the times of reporting and not the times 
of onset. 

As the disease began to spread, population mobility decreased. This was in part due to 
imposed regulations attempting to contain the disease but also due to the population’s fear 
of contagion. We reflect this in our model by decreasing the mobility value. Figure [3](b) 
shows the arrival times produced by our model when, as a strategy for slowing the spread, 
the mobility is reduced by 80%. Note that this reduction delays the arrival of EVD 
in each county by only a few weeks. This suggests that reducing the mobility of the 
individuals between counties will not stop the spread but will slow it sufficiently that 
other strategies can be developed and applied. Reducing mobility is also insufficient when 
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considering international transmission of the disease [lQj and more aggressive interventions 
are needed. We believe that an increase in both the percentage of infected individuals 
receiving hospitalisation in ETCs and the percentage of burials that follow procedures 
that do not transmit the disease are essential in containing the epidemic. 


Interventions and time to extinction 


To contain the disease and reduce its transmission we reduce mobility by 80%, increase 
the number of burials following procedures that do not transmit the disease, and increase 
the rate of hospitalisation in ETCs. Because health workers in ETCs have specialised 
training, we assume that the probability that they will be infected is greatly reduced and 
that the transmission coefficient /3# is decreased. A sufficiently rapid response to the 
EVD by the ETCs requires that (3h be decreased exponentially to a final value of 10” 3 , 
and hospitalisation 6 must be increased exponentially to reach 9—1. On the other hand, 
when changing local burial customs we assume that /3f decreases linearly and approaches 
zero. Changing local burial customs involves recruiting and training burial teams, takes 
a longer period of time, and is less aggressive than other kinds of intervention. This 
approach allows us to estimate an upper limit for the end of the epidemic, because we did 
not take into account other measures applied, such as, contact tracing, which could make 
the estimated end to the epidemic occur earlier. 

We apply these changes to simulate a two-month period, and the final result is that R 0 
decreases from 2.11 to i? 0 = 0.69, which is below the epidemic threshold. We consider two 
scenarios, (i) implementing the strategy beginning August 15 (the middle of the month 


indicated by WHO 


32] for the outbreak of EVD in West Africa) in which all symptomatic 


individuals are admitted to ETCs and safe burial procedures begin to apply, or (ii) im¬ 
plementing the same strategy, but beginning July 15 in order to study how delaying the 
implementation of the strategies affected containment. Our goal is to demonstrate that 
if the international response had been more rapid, the spreading disease would have been 
contained with a 50% probability by early March 2015 instead of the end of May 2015. 
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FIG. 4: Evolution of the number of cases (black) and deaths (red) when a reduction of 80% 
in the mobility rates is applied. The value of fin decreases exponentially to reach the value 
10~ 3 and (Up decreases linearly to reach a 0% of their original values. Also the hospitalisation 
fraction increases exponentially to reach 8 = 1. All reductions in the transmission coefficient 
were applied during two months, for (a) beginning at July 15th and (b) August 15th. Solid 
lines were obtained from the evolution equations (fil l 10 11 and the symbols are the data. The 
box plots show the median, the 25th and 75th percentile and 95% confidence interval of the 
median, obtained from the stochastic simulations. Figures c) and d) are the distribution of 
time to extinction of the EVD epidemic obtained from the stochastic simulations, when the 
strategy is applied from middle July and from middle August 2014, respectively. We show these 
distributions from December 1st, 2014 to December 1st, 2015. 


Figures 01(a) and 01(b) show that reducing the number of cases produced in hospitals 
and funerals reduces the cumulative number of cases to a plateau lower than the one 


predicted when no strategies are applied 


331 ]. Figure 01(a) shows that if our strategy had 


been applied in the middle of July the cumulative number of cases and deaths would 
have been approximately 80% lower than the reported number that resulted when the 
strategies were instituted in the middle of August. Figure 0)(b) shows that when we apply 
the strategy of our model to the actual mid-August starting time, it predicts (between 
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the 95% confidence interval of the median) the actual trend of cases and deaths reported 
in the WHO data in mid-March 2015. 

Our stochastic model allows us to quantify how the two different strategy implemen¬ 
tation times affect the extinction time of the EVD epidemic. Figures 2(c) and 2(d) show 
the extinction time distributions, i.e., when E = I = H = F = 0, when the strategy is 
implemented in July 2014 and August 2014, respectively (for the initial conditions, we 
use the cases provided by WHO for these dates). We find that the median of this distri¬ 
bution when the strategy is implemented in July is 6 March 2015 (with a 95% confidence 
interval from 5 January 2015 to 1 July 2015) and when it is implemented in August is 25 
May 2015 (with a 95% confidence interval from 28 March 2015 to 20 September 2015). 
Implementation in mid-August generated 8,000 cases of the disease, but an implemen¬ 
tation in mid-July would have reduced the time to disease extinction by three months 
and generated only 1,700 cases. The mid-August implementation faced a larger number 
of cases, the disease progression had a greater inertia against the strategy, and the cu¬ 
mulative number of cases required a longer time to go from an exponential regime to a 
subexponential regime. Thus if the health authorities and the international community 
had acted sooner the number of infected people would have been much lower. 


DISCUSSION 

In this manuscript we study the spreading of the Ebola virus using stochastic and 
deterministic compartmental models that incorporate the mobility of individuals between 
the counties in Liberia. We find that our model describes well the arrival of the disease into 
each of the counties, that reducing population mobility has little effect on geographical 
containment of the disease, and that reducing population mobility must be accompanied 
by other intervention strategies. We thus examine the effect of an intervention strategy 
that focuses on both an increase in safer hospitalisation and an increase in safer burial 
practices. Our study indicates that the intervention implemented in August 2014 reduced 
the total number of infected individuals significantly when compared to a scenario in 
which there is no strategy implementation, and it predicts that the epidemic will be 
extinct by mid-spring 2015. We also use our model to consider the difference in outcome 
had the strategy been implemented one month earlier. We find that the cumulative 
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number of cases and deaths would have been significantly lower and that the epidemic 
would have ended three months earlier. This indicates that a rapid and early intervention 
that increases the hospitalisation and reduces the disease transmission in hospitals and 
at funerals is the most important response to any possible re-emerging Ebola epidemic. 

Although our model simplifies the dynamics of epidemic spreading, it provides an 
adequate picture of the evolution in the number of cases and deaths. In future research 
we will incorporate more aspects of population mobility and intervention strategies carried 
out by health authorities. This will enable us to describe in greater detail the evolution 
of an epidemic and the efficacy of different strategies. 

Finally, the methods used in this manuscript to study Liberia can also be applied 
to Guinea and Sierra Leone as soon as high quality epidemic data from those countries 
become available. Future work should include both countries in order to quantify the 
cases spreading from them into Liberia. 


METHODS 

Stochastic model 

We generate a stochastic compartmental model based on the Gillespie algorithm. At 
each iteration of the simulation we draw a random number r (which represents the waiting 
time until the next transition) from an exponential distribution with parameter A given 
by parameter 

Nt Nco N co Nco 

A = EEWEE(£ j + ^>W% (12) 

i =1 j =1 i =1 j =1 

Here the first term A j is the rate of transition between states i in county j given in 
Table SI, and the second term corresponds to the mobility rates given in Eq. (fTTT) with 
x = E and x = S. 

Calibration with the deterministic equations 

To estimate the transmission coefficients /h, /3n, and /3f we calibrate a system of 
differential equations using least-squares optimisation with the data of the total cases 
from Liberia in the March-August period Q], and we apply a temporal shift, which 
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we will explain below. We compute the least-square values using a set of parameters 
generated using Latin hypercube sampling (LHS) in the parameter space [0, l] 3 , which 
we divide into 10 6 cubes of the same size. For each cube we choose a random point as a 
candidate for (/?/, /3h, Pf) in order to compute the standard deviation between the data 
and the system of differential equations obtained from this point. At the beginning of 
the epidemic there are very few cases (infected individuals), thus the evolution of the 
disease is in a stochastic regime in which the dispersion of the number of new cases is 
comparable to its mean value (see Refs. 344361]). When the number of infected individuals 
increases to a certain level, however, the epidemic evolves toward a quasi-deterministic 
regime and the evolution of the states of the stochastic simulation is the same as the 
states obtained using the solution of the evolution equations (Eqs. HlflUl). Nevertheless, 
due to fluctuations in the initial stochastic regime, a random temporal displacement of 
the quasi-deterministic growth of the number of accumulated cases is generated. Thus to 
remove this stochastic temporal shift and to compare the three aspects—the simulations, 
the numerical solution, and the data—we set the initial time at t — 0 when the total 
number of cases is above a cutoff * 0Q. For tire calibration of tire transmission 
coefficients, we use s c = 200 (which corresponds to the cumulative number of cases after 
21 July, according to the WHO data [l|), and using the least square method we give 
the data above this cutoff 50% of the weight because we are assuming that above s c 
the evolution of the disease spreading is quasi-deterministic. Finally, after we compute 
the sum of square residuals for each point in the parameter space, we apply the Akaike 
information criterion (AIC) and average those candidates of (/?/, /3h, Pf) with a AIC 
difference A < 2 [Q to obtain a model-averaged estimate of the transmission coefficients. 
An alternative method for estimating the transmission coefficients using an exponential 
fitting is discussed in Supplementary Information: Calibration. We find that this fitting 
generates the same set of values of transmission coefficients than the method with a 
temporal shift s c . Additionally, in Supplementary Information: Sensitivity Analysis we 
discuss the sensitivity of the estimated values of the transmission coefficients when 6 and 
5 change. 

The mobility data for the Wesolowski model were provided by Flowminder 
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26 


38] 


and the total cumulative case data used to calibrate the model were those supplied in 
reports generated by WHO jfjj]. Note that in this work we do not calibrate the model to 
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cases in each comity because it was shown by Chowell et al. 39] that globally the number 
of cases grows exponentially but locally can be better approximated by a polynomial 
than by an exponential growth. This cannot be addressed using our model because 
mathematically a differential equation with constant rates only reproduces an exponential 
growth. 


Weitz and Dushoff 


40] recently demonstrated that calibration causes an identification 


problem, i.e., that many combinations of the coefficient transmission values reproduce 
the real evolution of the number of cases, which is compatible with our finding that 
the calibrated transmission coefficients are in a plane (see Supplementary Information: 
Calibration). This point should be addressed in future research. 


Estimation of R n 


In order to compute t 


and Diekmann et al. 


re reproduction number Rq , following van den Driessche et al. 28] 
29], we construct a next-generation matrix. 

First, using the Jacobian matrix of the system of Eqs. (1-10) we construct the “trans¬ 
mission matrix” F. and the “transition matrix” V , obtaining 
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with i 7^ j. 

Note that the mobility rates are only in the transition matrix. Using these matrices we 
construct the next generation matrix, defined as FV _1 . Finally, the reproduction number 
is given by the spectral ratio p of the next generation matrix, R 0 = /^(FV” 1 ), i.e., its 
highest eigenvalue. Note that when the mobility rates go to zero, Rq decreases, i.e., in 
this limit an infected individual in a given county cannot interact with people from other 
comities and can only transmit the disease to susceptible individuals in the same county. 
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Transition Parameters 

Value 

References 

Mean duration of the incubation period (1/a) 

7 days 

[41-43] 

Mean time from the onset to the hospitalisation (I /7 h) 

5 days 

M 

Mean duration from onset to death (1/^d) 

9.6 days 

m 

Mean time from onset to the end for the cured (I/ 7 /) 

10 days 

[43, 45] 

Mean time from death to traditional burial ( 1 / 7 ^) 

2 days 

M 

Proportion of cases hospitalised (6) 

50% 

[15] 

Fatality Ratio (<5) 

50 % 

[15] 

Mean time from hospitalization to end for cured (I/ 77 ) 

5 days 

M 

Mean time from hospitalization to dead (I /7 hd) 

4.6 days 

M 


TABLE I: Transition parameters used to calculate the transition rates in our epidemic 
model. Table describing the different parameters used to calculate the transition rates among 
the ten different compartmental states in our model. 
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SUPPLEMENTARY INFORMATION 


Calibration 

In the main text, using the least square method and a temporal shift for the number of 
cases (s c = 200) we obtained the coefficient parameters /3j, /3h and /3f- We can also obtain 
these transmission coefficients by measuring the parameter 77 that best fits an exponential 
growth I(t) ~ exp(r/t) in the deterministic regime. For the epidemic spreading in Liberia, 
we find that fitting with an exponential function on the number of cases from July 21st 
to August 15th, r] = 0.053 ± 0.003. To find the relation between 77 and the transmission 
parameters /?/, /3h, and /3p, we write the equation of the eigenvalue of the Jacobian J of 
the system of Eqs. (1)-(10) for the N co = 15 counties 

det^J^^H^F) -rj IdJ = 0 , (SI) 

where det is the determinant function, rj is the eigenvalue obtained from the fitting, Id 
is the identity matrix with 10JV CO rows and columns, where the factor 10 corresponds 
to the number of evolution equations for each county. Note that J is a function of the 
transmission coefficients. 

Equation (ISID sets the relationship between the exponential growth rate ( 77 ) and the 
transmission parameters, which is not linear when the mobility is taken into account. 
However, we can approximate this equation by neglecting the flow of individuals. This 
is the case because although mobility spreads the EVD throughout the country in the 
stochastic stage, when the disease reaches the deterministic regime its spreading is pri¬ 
marily due to infected individuals within each county and imported cases are no longer 
a relevant factor. With this approximation, and using 77 = 0.054, we obtain an equation 
of a plane (see Fig. ISTf that coincides with the triads of transmission coefficients that 
were obtained by using the least square fitting with a shift s c . Thus this method and the 
exponential fitting generate the same set of transmission parameters that reproduces the 
epidemic growth. 


18 


0.0 



FIG. SI: Values of (Pi, Ph, Pf) that fits the data obtained from i) the shift and least square 
method (points) and ii) fitting the exponential growth of number of cases (plane). 

Note that in the main text we use s c = 200 to calibrate the data. On the other hand, if 
we use s c = 100 (near the 2 July date) we obtain P / = 0.09 (0, 0.23), Ph = 0.31 (0, 0.78), 
and Pf = 0.46 (0,0.99), and if we use s c = 300 (near the 28 July date) we obtain 
Pr = 0.11 (0,0.27), p H = 0.38 (0,0.94) and p F = 0.46 (0,0.99). 
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Table of the transitions 


Transition 

Transition rate (Aj) 

(S,E) (S — 1, E + 1) 

MPf S I + pH s H + p F S F) 

(E, Idh ) —> (E — 1, Idh + 1) 

a 8 5E 

(E, Idnh ) —> (E - 1, Idnh + 1) 

a (1-8) 5E 

(E, Idh) —> (E — 1, Irh + 1) 

a 8 (l — S) E 

{E, Irnh) (E — 1, Irnh + 1) 

a (1-8) (1-8) E 

(Idh, H d ) -t ( Idh — 1, H d + 1) 

7 h Idh 

(Idnh, E) —>■ (Idnh - 1, F + 1) 

7 d Idnh 

(Irh, Hr) —> (I RH - 1, Hr + 1) 

7 h Irh 

(Irnh, R) —> (Irnh - 1, R + 1) 

7 1 Irnh 

(H d ,F)^(H d -1,F+1) 

Ihd Hd 

(Hr, R) -»• (H r -1,R+1) 

1HI Hr 

(F,R)->(F-1,R+1) 

7 F F 


TABLE SI: Table of the transition with their respective transition rates for our 
model. Table representing the transition rates between different compartmental states in our 
model. The capital letters represents number of: susceptible individuals (S), number of exposed 
individuals (E), individuals infected who will be hospitalised and die ( Idh ), individuals infected 
who won’t be non hospitalised and will die ( Idnh )> individuals infected who will be hospitalised 
and recovered ( Irh ), individuals infected who won’t be hospitalised and will recover ( Irnh), 
individuals hospitalised who will die (Hd ), individuals hospitalised who will recover (Hu). Here 
R is the number of individuals cured or dead and F is the number of individuals in the funerals 
who will have unsafe burials and can infect. Here Pi, Ph and Pf are the transmission coefficients 
in the community in the hospital and in the funerals respectively, 5 is the fatality ratio and 8 
the fraction of the hospitalised ones. The inverse of the mean time period of the incubation is 
1/a. The mean time period from symptoms to hospitalisation is I/ 7 h, from symptoms for non 
hospitalised individuals to dead is I/'Jd, from symptoms for hospitalised individuals to dead 
is 1 /^hD, from symptoms for hospitalised individuals to recovery is 1 /jhi and from dead to 
recover is I/ 7 p- The flow of mobility for individuals in county i —>• j is explained in Eq. (11). 
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Sensitivity analysis 


Because the epidemic evolution in this model is affected by many parameters, small 
changes in their values could significantly impact the model’s output. Here we analyse 
how changing the hospitalisation ratio 9 and the death ratio 8 affects the estimation of /3j, 
(3h, Pf , and Ro , and also affects the evolution of the cumulative cases when the strategy 
is implemented in August as described in the main text. 

For 9 and 5 = 0.40, 0.50, 0.60 we show the values of P I} p Hl ftp, and R 0 , using the 
least square method and the Akaike average, as explained in the Methods section. 


e 

Pi 

Ph 

Pf 

Ro 

0.40 

0.13 

0.37 

0.40 

2.24 (1.99,2.28) 

0.50 

0.14 

0.29 

0.40 

2.11 (1.88,2.71) 

0.60 

0.14 

0.25 

0.39 

2.05 (1.92,2.28) 


TABLE S2: Values of the transmission coefficients and Ro for 0 = 0.40, 0.50, 0.60. Here 5 = 0.5 


5 

Pi 

Ph 

Pf 

Ro 

0.40 

0.14 

0.30 

0.42 

2.18 (1.97,2.16) 

0.50 

0.14 

0.29 

0.40 

2.11 (1.88,2.71) 

0.60 

0.13 

0.29 

0.38 

2.05 (1.96,2.63) 


TABLE S3: Values of the transmission coefficients and R$ for 5 = 0.40, 0.50, 0.60. Here 6 = 0.5 


Table [S2] shows that as 9 increases 50% from 9 = 0.40 to 9 = 0.60, /3h decreases 
30%. In contrast, /3j and ftp remain almost constant and the reproductive number does 
not change significantly. On the other hand, Table [S3] shows that when 8 increases from 
6 = 0.40 to 8 = 0.60, the transmission coefficients and Rq change less than 10%, however 
it remains inside the interval. Thus this model is more sensitive to changes in 9 than in 
6. Figure [S2] plots the number of cumulative cases as a function of time obtained from 
the stochastic model, compares the results with WHO data, and shows good agreement 
between the simulations and the real data. 
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FIG. S2: Evolution of the cumulative number of cases in Liberia with 100 realisations (gray lines) 
and the data (symbols) with a temporal shift using s c = 200. Figures (a) and (b) correspond 
to 0 = 0.40 and 6 = 0.60 respectively, with 5 = 0.5. The transmission coefficients used, were 
obtained from table [S2l Figures (c) and (d) correspond to 5 = 0.40 and 6 = 0.60, respectively, 
with 6 = 0.50. The transmission coefficients used, were obtained from table IS3l 


To evaluate how changing the parameter values alters the effectiveness of the inter¬ 
vention strategy as explained in the main text, Fig. [S3] plots the cumulative number of 
infected individuals when the strategy is implemented in mid-August for different values 
of 0 and 5. 
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FIG. S3: Evolution of the cumulative number of cases (black) and deaths (red) when it is applied 
the strategy from August, as it was explained in the main text, for different values of 9 and 6. 
In the figure (a), 8 = 0.50 and 8 = 0.50 (solid line), 9 = 0.40 (dotted line and filled box plots) 
and 9 = 0.60 (dashed line and open box plots). In figure (b) 9 = 0.50 and 6 = 0.50 (solid line), 
S = 0.40 (dotted line and filled box plots) and 9 = 0.60 (dashed line and open box plots). All 
the curves were obtained by integrating the evolution equations (1)-(10), and the box plots were 
obtained from the stochastic simulations. 

Figure [S3] shows that although the evolution of the number of cases is not sensitive to 
variations in 9 and S, the evolution of number of deaths is sensitive to variations in 6 (see 
FigJ53b). This is the case because the final number of deaths is proportional to 6, and in 
our strategy /3p changes but 5 remains constant. 


[1] The New York Times / Sack, K., Fink, S., Belluck, P., Nossiter, 
A. & Berehulak, D. How Ebola Roared Back. (2014). Available at: 
http://www.nytimes.com/2014/12/30/health/how-ebola-roared-back.html, (Accessed: 
4th February 2015). 

[2] Frieden, T. R., Damon, I., Bell, B. P., Kenyon, T. Sz Nichol, S. Ebola 2014 - New Challenges, 
New Global Response and Responsibility. N. Engl. J. Med. 371, 1177-1180 (2014). 

[3] World Health Organization. Ebola response roadmap - Situation report (2015). Available 
at: http://www.who.int/csr/disease/ebola/situation-reports/en/ (Accessed: 30th March 
2015). 


23 













[4] Center for Disease Control and Prevention. 2014 Ebola Outbreak in West Africa (2014). 

Available at: http://www.cdc.gov/vhf/ebola/outbreaks/2014-west-africa/index.html, (Ac¬ 
cessed: 1st January 2015). 

[5] Aljazeera America / A. Mukpo. The biggest concern of the Ebola outbreak is political, not 
medical (2014). Available at: http://america.aljazeera.com/opinions/2014/8/ebola-virus- 
liberiasierraleonepolitics.html, (Accessed: 4th February 2015). 

[6] The Guardian / A. Konneh. Ebola isn’t just a health crisis - it’s a social and economic one 
too. (2014). Available at: http://www.theguardian.com/commentisfree/2014/oct/10/ebola- 
liberia-catastrophe-generation-poverty, (Accessed: 4th February 2015). 

[7] Global Issues / T. Brewer. Water and Sanitation Report Card: Slow Progress, Inadequate 
Funding. (2014). Available at: http://www.globalissues.org/news/2014/ll/24/20342, (Ac¬ 
cessed: 4th February 2015). 

[8] The Washington Post / L.H. Sun, B. Dennis and J. Achenbach. Ebola’s lessons, 

painfully learned at great cost in dollars and human lives (2014). Available at: 
http://www.washingtonpost.com/national/health-science/ebolas-lessons-painfully-learned-at-great-cos 
(Accessed: 4th February 2015). 

[9] World Bulletin / News Desk. 2014: Liberia’s Ebola nightmare (2014). Available at: 
http://www.worldbulletin.net/news/151832/2014-liberias-ebola-nightmare, (Accessed: 4th 
February 2015). 

[10] Scientific American / S. Yasmin. Ebola Infections Fewer Than Predicted by Disease Mod¬ 
els (2014). Available at: http://www.scientificamerican.com/article/ebola-infections-fewer- 
than-predicted-by-disease-models/, (Accessed: 4th February 2015). 

[11] News 24. Ebola: Bad data an issue (2014). Available at: 

http://www.news24.com/Africa/News/Ebola-Bad-data-an-issue-20141208, (Accessed: 

4th February 2015). 

[12] The New Zealand Herald. In Ebola outbreak, bad data adds another problem (2014). Avail¬ 
able at: http://www.nzherald.co.nz/world/news/article.cfm7c_id=2&objectid=11371063, 

(Accessed: 4th February 2015). 

[13] World Health Organization. Ebola Virus Disease Outbreak Response Plan in West Africa 
(2014). Available at: http://www.who.int/csr/disease/ebola/outbreak-response-plan/en/, 

(Accessed: 4th February 2015). 


24 


[14] Legrand, J., Grais, R., Boelle, P., Valleron, A. & Flahault, A. Understanding the dynamics 
of Ebola epidemics. Epidemiol. Infect. 135, 610-621 (2007). 

[15] Rivers, C. M., Lofgren, E. T., Marathe, M., Eubank, S. & Lewis, B. L. Modeling the 
Impact of Interventions on an Epidemic of Ebola in Sierra Leone and Liberia. PL OS Curr.: 
Outbreaks. DOI: 10.1371/currents.outbreaks.fd38dd85078565450b0be3fcd78f5ccf (2014). 

[16] Gomes, M. F. et al. Assessing the international spreading risk associated with the 
2014 West African Ebola outbreak. PLOS Curr.: Outbreaks. DOI: 10.1371/cur¬ 
rents.outbreaks.Cd818f63d40e24aef769dda7df9e0da5 (2014). 

[17] Balcan D. et al. Seasonal transmission potential and activity peaks of the new influenza 
A(H1N1): a Monte Carlo likelihood analysis based on human mobility. BMC Med. 7, 45 
(2009). 

[18] Balcan D. et al. Multiscale mobility networks and the spatial spreading of infectious dis¬ 
eases. Proc. Natl. Acad. Sci. U.S.A. 106 , 21484-9 (2009). 

[19] Poletto, C. et al. Assessing the impact of travel restrictions on international spread of the 
2014 West African Ebola epidemic. Euro. Surveill. 19, 20936 (2014). 

[20] Merler, S. et al. Spatiotemporal spread of the 2014 outbreak of Ebola virus disease in 
Liberia and the effectiveness of non-pharmaceutical interventions: a computational mod¬ 
elling analysis. Lancet Infect. Dis. 15, 204-211 (2015). 

[21] Wesolowski, A. et al. Commentary: Containing the Ebola outbreak-the potential 
and challenge of mobile network data. PLOS Curr.: Outbreaks. DOI: 10.1371/cur¬ 
rents.outbreaks.0177e7fcf52217b8b634376e2f3efc5e (2014). 

[22] Liberia Institute of Statistics and Geo-Information Services. (2009). 2008 Population and 
Housing Census Final Results [Online]. Monrovia, Liberia: Government of the Republic of 
Liberia. (2009). Available at: http://www.lisgis.net/pgjing/NPHC2008 Final Report.pdf, 
(Accessed: 23rd February 2015). 

[23] Gonzalez, M. C., Hidalgo, C. A. & Barabasi, A.-L. Understanding individual human mo¬ 
bility patterns. Nature 453, 779-782 (2008). 

[24] Tatem, A. J. et al. The use of mobile phone data for the estimation of the travel patterns 
and imported Plasmodium falciparum rates among Zanzibar residents. Malar. J. 8, 287 
(2009). 

[25] Wesolowski, A. et al. Quantifying travel behavior for infectious disease research: a compar- 


25 



ison of data from surveys and mobile phones. Sci. Rep. 4, 5678; DOI: 10.1038/srep05678 
(2014). 

[26] Flowminder Foundation ( 2014). Available at: http://www.flowminder.org/, (Accessed: 4th 
February 2015). 

[27] Colizza, V. & Vespignani, A. Invasion threshold in heterogeneous metapopulation networks. 
Phys. Rev. Lett. 99 , 148701 (2007). 

[28] Van den Driessche, P. & Watmough, J. Reproduction numbers and sub-threshold endemic 
equilibria for compartmental models of disease transmission. Math. Biosci. 180 , 29-48 
( 2002 ). 

[29] Diekmann, O., Heesterbeek, J. & Roberts, M. The construction of next-generation matrices 
for compartmental epidemic models. J. R. Soc. Interface 7 , 873-885 (2010). 

[30] Anderson, R. M. & May, R. M. Infectious Diseases of Humans: Dynamics and Control 
(Oxford University Press, Oxford, 1992). 

[31] Daily Observer. 2 Of 5 Test Positive For Ebola In Liberia. (2014). Available 
at: http://www.liberianobserver.com/security-health/2-5-test-positive-ebola-liberia, (Ac¬ 
cessed: 4th February 2015). 

[32] WHO Statement on the Meeting of the International Health Regulations Emergency 

Committee regarding the 2014 Ebola outbreak in West Africa.( 2014). Available at: 

http://www.who.int/mediacentre/news/statements/2014/ebola-20140808/en/, (Accessed: 
4th February 2015). 

[33] Meltzer, M. I. et al. Estimating the future number of cases in the ebola epidemic—Liberia 
and Sierra Leone, 2014-2015. MMWR Surveill Summ 63, 1-14 (2014). 

[34] Volz, E. SIR dynamics in random networks with heterogeneous connectivity. J. Math. Biol. 
56, 293-310 (2008). 

[35] Valdez, L. D., Macri, P. A. & Braunstein, L. A. Temporal Percolation of the Susceptible Net¬ 
work in an Epidemic Spreading. PLoS ONE 7 , e44188; DOI: 10.1371/journal.pone.0044188 
( 2012 ). 

[36] Barthelemy, M., Barrat, A., Pastor-Satorras, R. & Vespignani, A. Dynamical patterns of 
epidemic outbreaks in complex heterogeneous networks. J. Theor. Biol. 235 , 275 - 288 
(2005). 

[37] Burnham, K. P. & Anderson, D. R. Model selection and multimodel inference: a practical 


26 


information-theoretic approach (Springer, 2002). 

[38] The WorldPop project. Mobile phone flow maps (2014). Available at: 
http://www.worldpop.org.uk/ebola/, (Accessed: 4th February 2015). 

[39] Chowell, G., Viboud, C., Hyman, J. M. & Simonsen, L. The Western Africa Ebola virus 
disease epidemic exhibits both global exponential and local polynomial growth rates. PLOS 
Curr.: Outbreaks. DOI: 10.1371/currents.outbreaks.8b55f4bad99ac5c5db3663e916803261 
(2015). 

[40] Weitz, J. S. & Dushoff, J. Modeling Post-death Transmission of Ebola: Challenges for 
Inference and Opportunities for Control. Sci. Rep. 5, 8751; DOI: 10.1038/srep08751 (2015). 

[41] Bwaka, M. A. et al. Ebola hemorrhagic fever in Kikwit, Democratic Republic of the Congo: 
clinical observations in 103 patients. J. Infect. Dis. 179, S1-S7 (1999). 

[42] Ndarnbi, R. et al. Epidemiologic and clinical aspects of the Ebola virus epidemic in Mosango, 
Democratic Republic of the Congo, 1995. J. Infect. Dis. 179, S8-S10 (1999). 

[43] Dowell, S. F. et al. Transmission of Ebola hemorrhagic fever: a study of risk factors in 
family members, Kikwit, Democratic Republic of the Congo, 1995. J. Infect. Dis. 179, 
S87-S91 (1999). 

[44] Khan, A. S. et al. The reemergence of Ebola hemorrhagic fever, Democratic Republic of 
the Congo, 1995. J. Infect. Dis. 179, S76-S86 (1999). 

[45] Rowe, A. K. et al. Clinical, virologic, and immunologic follow-up of convalescent Ebola 
hemorrhagic fever patients and their household contacts, Kikwit, Democratic Republic of 
the Congo. J. Infect. Dis. 179, S28-S35 (1999). 

Acknowledgements 

HES thanks the NSF (grants CMMI 1125290 and CHE-1213217) and the Keck Foundation 
for financial support. LDV and LAB wish to thank to UNMdP and FONCyT (Piet 
0429/2013) for financial support. 

Contributions 

LDV, HHAR, HES and LAB designed the research, analyzed data, discussed results, 
and contributed to writing the manuscript. LDV implemented and performed numerical 
experiments and simulations. 

Additional information 


27 


Competing financial interests: The authors declare no competing financial interests. 


